library(ggplot2)
library(haven)
library(list)
library(dplyr)
library(stargazer)
library(estimatr)
library(ggpubr)
library(scales)
library(xtable)
library(MASS)
library(ordinal)


#####################################################################################################################
# Replication code for "Corrupted Estimates? Response Bias in Citizen Surveys on Corruption", Political Behavior
# Author: Mattias Agerberg, University of Gothenburg, mattias.agerberg@gu.se
# SUPPLEMENTARY MATERIALS
#####################################################################################################################


#####################################################################################################################
# Paper wd
setwd("") 

# Data
data <- read.csv("data_RO.csv") 

#####################################################################################################################
# ADDITIONAL DATA MANAGEMENT + TESTS AND DESCRIPTIVES
#####################################################################################################################
attach(data)

# Direct bribe question
data$bribe[Bribe_direct=="Nu"] <- 0
data$bribe[Bribe_direct=="Prefer sa nu raspund"] <- 0
data$bribe[Bribe_direct=="Da"] <- 1

# Reverse code Corruption change and Economy - high values: corruption increase, worse economy
data$corr_soc <- data$corr_soc*(-1)+6
data$econ <- data$econ*(-1)+6

# Ranked response time variable
data$rank_response <- ntile(data$Duration__in_seconds_,100)

# List variables
list <- cbind(List_nonsensitive,List_sensitive)
data$items[is.na(list[,1]==F) | is.na(list[,2]==F)] <- rowSums(list, na.rm = TRUE)
data$items <- as.integer(data$items)

data$list_treatment <- NA
data$list_treatment[is.na(List_nonsensitive)==F] <- 0
data$list_treatment[is.na(List_sensitive)==F] <- 1
data$list_treatment <- as.integer(data$list_treatment)

data_complete <- data[,c("items","list_treatment","gov_support_01","female","city_over200k","university_ed","Age","income_top20","rank_response")]
data_complete <- data_complete[complete.cases(data_complete),]

# TABLE 9
# Test for design effect + table - tablenotes were added using "threeparttable" in LaTeX
test.value.list <- ict.test(data$items, data$list_treatment, J = 4, gms = TRUE, pi.table = T)
print(test.value.list)

stargazer(test.value.list$pi.table, title = "Respondent types, estimated proportions", notes = "The table shows the estimated proportion of respondent types characterized by the total number of affirmative answers to the control questions, $Y$, and the truthful answer for the sensitive item $Z$ (1 indicates affirmative and 0 represents negative). Standard errors are also provided for each estimated proportion.")

# TABLE 3
# Descriptives (sample)
descriptives <- data[,c("female","city_over200k","university_ed","Age","Income","Household")]
descriptives$Household <- as.numeric(descriptives$Household)
descriptives$Income[data$Income>1000000] <- NA
median(descriptives$Income, na.rm = TRUE)
IQR(descriptives$Income, na.rm = TRUE)
stargazer(descriptives, summary.stat = c("mean","sd","min","max","n"), 
          covariate.labels = c("Female respondent","City with over 200,000 inhabitants",
                               "University education","Age","Household income (Lei per month)","Persons in household"), digits = 2,
          title = "Sample characteristics", 
          notes = "Note: Some extreme (probably miscoded) outliers were excluded from the `Income' variable. For the Income variable, the median household income is reported, and the interquartile range (IQR) is reported in the `St. Dev.' column.", notes.append = TRUE)

# TABLE 4
# Descriptives (attitudes)
descriptives2 <- data[,c("gov_supp","gov_support_01","pol_interest","corr_soc","econ","corr_pol","corr_worry","bribe")]
stargazer(descriptives2, summary.stat = c("mean","sd","min","max","n"), 
          covariate.labels = c("Government attitude (5=favoring)","Government supporter (attitude+party)", "Political interest (4=high)",
                               "Corruption increase","Worse economy","Corruption in politics",
                               "Corruption worry","Direct bribe question (1=yes)"), digits = 2,
          title = "Respondent attitudes")

# FIGURE 5
# Parties
data$party <- factor(data$party, labels = c("PSD","ALDE","UDMR","PNL","USR","PMP","Other party","Would not vote"))

ggplot(data, aes(x = factor(party))) +  
  geom_bar(aes(y = (..count..)/sum(..count..)), colour="black", fill="grey20", alpha = 0.75) +
  theme_minimal() +
  labs(x="", y="Share") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggsave("parties.pdf")

# Question order tests
data$control_groups <- factor(data$control_groups, 
                              labels = c("Corr. politics \n control","Corr. increase \n control", "Econ. worse \n control", "Corr. worry \n control"))

ordertest.1 <- glm(gov_party~control_groups, data = data, family = "binomial")
anova(ordertest.1, test = "Chisq")

data$gov_supp <- factor(data$gov_supp)
ordertest.02 <- polr(gov_supp~1, data = data, Hess = TRUE)
ordertest.2 <- polr(gov_supp~control_groups, data = data, Hess = TRUE)
anova(ordertest.02,ordertest.2, test = "Chisq")

#####################################################################################################################
# LIST ANALYSIS, SUPPLEMENTARY
#####################################################################################################################
# FIGURE 8
# List analysis, subgroups fit
fit.list <- ictreg(items ~ gov_support_01 + female + city_over200k + university_ed + income_top20 + Age,
                   data = data_complete, treat = "list_treatment", J = 4, method = "nls")

fit.sens <- glm(bribe ~ gov_support_01 + female + city_over200k + university_ed + income_top20 + Age,
                data = data, family = binomial("logit"))

# List analysis, Education
data_complete_nouni <- data_complete_uni <- data_complete
data_complete_uni[, "university_ed"] <- 1
data_complete_nouni[, "university_ed"] <- 0

avg.pred.uni <- predict(fit.list,
                        newdata = data_complete_uni, avg = TRUE, interval = "confidence")

avg.pred.nouni <- predict(fit.list,
                          newdata = data_complete_nouni, avg = TRUE, interval = "confidence")

newdata2 <- with(data, data.frame(female = mean(female,na.rm = TRUE), gov_support_01 = mean(gov_support_01), city_over200k = mean(city_over200k), 
                                  university_ed = 1, income_top20 = mean(income_top20), Age = mean(Age)))
pred.sens.uni <- predict(fit.sens, newdata = newdata2, type = "response", se.fit = TRUE, avg = TRUE)

newdata2 <- with(data, data.frame(female = mean(female,na.rm = TRUE), gov_support_01 = mean(gov_support_01), city_over200k = mean(city_over200k), 
                                  university_ed = 0, income_top20 = mean(income_top20), Age = mean(Age)))
pred.sens.nouni <- predict(fit.sens, newdata = newdata2, type = "response", se.fit = TRUE, avg = TRUE)

# Plot by education
pred.mf <- matrix(nrow = 4,ncol = 5)
pred.mf[c(1,3),4] <- "University degree"
pred.mf[c(2,4),4] <- "No uni. degree"
pred.mf[1:2,5] <- "List"
pred.mf[3:4,5] <- "Direct"
for (i in 1:3) {
  pred.mf[1,i] <- as.numeric(avg.pred.uni$fit[[i]])
}
for (i in 1:3) {
  pred.mf[2,i] <- as.numeric(avg.pred.nouni$fit[[i]])
}
pred.mf[3,1] <- pred.sens.uni[[1]]
pred.mf[3,2] <- pred.sens.uni[[1]]-1.96*pred.sens.uni[[2]]
pred.mf[3,3] <- pred.sens.uni[[1]]+1.96*pred.sens.uni[[2]]
pred.mf[4,1] <- pred.sens.nouni[[1]]
pred.mf[4,2] <- pred.sens.nouni[[1]]-1.96*pred.sens.nouni[[2]]
pred.mf[4,3] <- pred.sens.nouni[[1]]+1.96*pred.sens.nouni[[2]]

pred.mf <- data.frame(pred.mf)
pred.mf$X1 <- as.numeric(levels(pred.mf$X1))[pred.mf$X1]
pred.mf$X2 <- as.numeric(levels(pred.mf$X2))[pred.mf$X2]
pred.mf$X3 <- as.numeric(levels(pred.mf$X3))[pred.mf$X3]

edu.list.gr <- ggplot(pred.mf,aes(x=X4,y=X1,colour=X5)) +
  geom_errorbar(aes(ymin=X2, ymax=X3), width=.1,position = position_dodge(width = 0.2), size = 0.25) +
  geom_point(position = position_dodge(width = 0.2),size=1.3) +
  scale_y_continuous(limits=c(0, 0.83)) +
  labs(y="",x="",colour="") +
  theme_minimal(base_size = 12) + scale_colour_grey(start=0.5, end=0.2) +
  theme(axis.text.x = element_text(size=6.5), axis.title.y = element_text(size=5.5, vjust = 1), 
        axis.text.y = element_text(size=5.3), legend.text = element_text(size=6.5),
        panel.grid.minor = element_line(size = 0.15), panel.grid.major = element_line(size = 0.3)) 

# List analysis, City
data_complete_nocity <- data_complete_city <- data_complete
data_complete_city[, "city_over200k"] <- 1
data_complete_nocity[, "city_over200k"] <- 0

avg.pred.city <- predict(fit.list,
                         newdata = data_complete_city, avg = TRUE, interval = "confidence")

avg.pred.nocity <- predict(fit.list,
                           newdata = data_complete_nocity, avg = TRUE, interval = "confidence")

newdata2 <- with(data, data.frame(female = mean(female,na.rm = TRUE), gov_support_01 = mean(gov_support_01), city_over200k = 1, 
                                  university_ed = mean(university_ed), income_top20 = mean(income_top20), Age = mean(Age)))
pred.sens.city <- predict(fit.sens, newdata = newdata2, type = "response", se.fit = TRUE, avg = TRUE)

newdata2 <- with(data, data.frame(female = mean(female,na.rm = TRUE), gov_support_01 = mean(gov_support_01), city_over200k = 0, 
                                  university_ed = mean(university_ed), income_top20 = mean(income_top20), Age = mean(Age)))
pred.sens.nocity <- predict(fit.sens, newdata = newdata2, type = "response", se.fit = TRUE, avg = TRUE)

# Plot by city/non-city
pred.mf <- matrix(nrow = 4,ncol = 5)
pred.mf[c(1,3),4] <- "City (over 200 000)"
pred.mf[c(2,4),4] <- "No city"
pred.mf[1:2,5] <- "List"
pred.mf[3:4,5] <- "Direct"
for (i in 1:3) {
  pred.mf[1,i] <- as.numeric(avg.pred.city$fit[[i]])
}
for (i in 1:3) {
  pred.mf[2,i] <- as.numeric(avg.pred.nocity$fit[[i]])
}
pred.mf[3,1] <- pred.sens.city[[1]]
pred.mf[3,2] <- pred.sens.city[[1]]-1.96*pred.sens.city[[2]]
pred.mf[3,3] <- pred.sens.city[[1]]+1.96*pred.sens.city[[2]]
pred.mf[4,1] <- pred.sens.nocity[[1]]
pred.mf[4,2] <- pred.sens.nocity[[1]]-1.96*pred.sens.nocity[[2]]
pred.mf[4,3] <- pred.sens.nocity[[1]]+1.96*pred.sens.nocity[[2]]

pred.mf <- data.frame(pred.mf)
pred.mf$X1 <- as.numeric(levels(pred.mf$X1))[pred.mf$X1]
pred.mf$X2 <- as.numeric(levels(pred.mf$X2))[pred.mf$X2]
pred.mf$X3 <- as.numeric(levels(pred.mf$X3))[pred.mf$X3]

city.list.gr <- ggplot(pred.mf,aes(x=X4,y=X1,colour=X5)) +
  geom_errorbar(aes(ymin=X2, ymax=X3), width=.1,position = position_dodge(width = 0.2), size = 0.25) +
  geom_point(position = position_dodge(width = 0.2),size=1.3) +
  scale_y_continuous(limits=c(0, 0.83)) +
  labs(y="",x="",colour="") +
  theme_minimal(base_size = 12) + scale_colour_grey(start=0.5, end=0.2) +
  theme(axis.text.x = element_text(size=6.5), axis.title.y = element_text(size=5.5, vjust = 1), 
        axis.text.y = element_text(size=5.3), legend.text = element_text(size=6.5),
        panel.grid.minor = element_line(size = 0.15), panel.grid.major = element_line(size = 0.3)) 

# List analysis, Age
data_complete_age50 <- data_complete_age25 <- data_complete
data_complete_age25[, "Age"] <- 25
data_complete_age50[, "Age"] <- 50

avg.pred.age25 <- predict(fit.list,
                          newdata = data_complete_age25, avg = TRUE, interval = "confidence")

avg.pred.age50 <- predict(fit.list,
                          newdata = data_complete_age50, avg = TRUE, interval = "confidence")

newdata2 <- with(data, data.frame(female = mean(female,na.rm = TRUE), gov_support_01 = mean(gov_support_01), city_over200k = mean(city_over200k), 
                                  university_ed = mean(university_ed), income_top20 = mean(income_top20), Age = 25))
pred.sens.age25 <- predict(fit.sens, newdata = newdata2, type = "response", se.fit = TRUE, avg = TRUE)

newdata2 <- with(data, data.frame(female = mean(female,na.rm = TRUE), gov_support_01 = mean(gov_support_01), city_over200k = mean(city_over200k), 
                                  university_ed = mean(university_ed), income_top20 = mean(income_top20), Age = 50))
pred.sens.age50 <- predict(fit.sens, newdata = newdata2, type = "response", se.fit = TRUE, avg = TRUE)

# Plot by young/old
pred.mf <- matrix(nrow = 4,ncol = 5)
pred.mf[c(1,3),4] <- "Age: 25"
pred.mf[c(2,4),4] <- "Age: 50"
pred.mf[1:2,5] <- "List"
pred.mf[3:4,5] <- "Direct"
for (i in 1:3) {
  pred.mf[1,i] <- as.numeric(avg.pred.age25$fit[[i]])
}
for (i in 1:3) {
  pred.mf[2,i] <- as.numeric(avg.pred.age50$fit[[i]])
}
pred.mf[3,1] <- pred.sens.age25[[1]]
pred.mf[3,2] <- pred.sens.age25[[1]]-1.96*pred.sens.age25[[2]]
pred.mf[3,3] <- pred.sens.age25[[1]]+1.96*pred.sens.age25[[2]]
pred.mf[4,1] <- pred.sens.age50[[1]]
pred.mf[4,2] <- pred.sens.age50[[1]]-1.96*pred.sens.age50[[2]]
pred.mf[4,3] <- pred.sens.age50[[1]]+1.96*pred.sens.age50[[2]]

pred.mf <- data.frame(pred.mf)
pred.mf$X1 <- as.numeric(levels(pred.mf$X1))[pred.mf$X1]
pred.mf$X2 <- as.numeric(levels(pred.mf$X2))[pred.mf$X2]
pred.mf$X3 <- as.numeric(levels(pred.mf$X3))[pred.mf$X3]

age.list.gr <- ggplot(pred.mf,aes(x=X4,y=X1,colour=X5)) +
  geom_errorbar(aes(ymin=X2, ymax=X3), width=.1,position = position_dodge(width = 0.2), size = 0.25) +
  geom_point(position = position_dodge(width = 0.2),size=1.3) +
  scale_y_continuous(limits=c(0, 0.83)) +
  labs(y="",x="",colour="") +
  theme_minimal(base_size = 12) + scale_colour_grey(start=0.5, end=0.2) +
  theme(axis.text.x = element_text(size=6.5), axis.title.y = element_text(size=5.5, vjust = 1), 
        axis.text.y = element_text(size=5.3), legend.text = element_text(size=6.5),
        panel.grid.minor = element_line(size = 0.15), panel.grid.major = element_line(size = 0.3)) 

# Combine list graphs, appendix
list.combined2 <- ggarrange(age.list.gr, city.list.gr, edu.list.gr,
                            ncol = 3, nrow = 1, common.legend = TRUE, legend = "bottom")
list.combined2 <- annotate_figure(list.combined2, left = text_grob("Estimated proportion",size = 7,rot = 90,hjust = 0.05,vjust = 3))
ggsave("listsubgroups2.pdf",width=18,height=8, units = "cm")

# FIGURE 9
# Difference-in-means, response time
fit.list <- ictreg(items ~ 1,
                   data = data[data$rank_response>5,], treat = "list_treatment", J = 4, method = "lm")

fit.sens <- glm(bribe ~ 1,
                data = data[data$rank_response>5,], family = binomial("logit"))

avg.pred.social.desirability <- predict(fit.list,
                                        direct.glm = fit.sens, se.fit = TRUE)

pred.bias <- as.matrix(avg.pred.social.desirability[[1]])
pred.bias <- data.frame(pred.bias)
pred.bias$est <- c("List \n estimate","Direct \n estimate","Difference")
pred.bias$est <- factor(pred.bias$est)
pred.bias$est = factor(pred.bias$est,levels(pred.bias$est)[c(2,3,1)])

pct5.gr <- ggplot(pred.bias,aes(x=est,y=fit,colour=est)) +
  geom_errorbar(aes(ymin=lwr, ymax=upr), width=.1,position = position_dodge(width = 0.2), size = 0.25) +
  geom_point(position = position_dodge(width = 0.2),size=1.3) +
  scale_y_continuous(limits=c(0,0.5)) +
  labs(y="",x="", title="5th pctl < X") +
  theme_minimal(base_size = 12) + scale_colour_grey(start=0.5, end=0.2) +
  theme(legend.position="none", axis.title.y = element_text(vjust = 2.4),axis.text.x = element_text(size=6.5),
        axis.text.y = element_text(size=5.3), plot.title = element_text(size=9))

fit.list <- ictreg(items ~ 1,
                   data = data[data$rank_response>10,], treat = "list_treatment", J = 4, method = "lm")

fit.sens <- glm(bribe ~ 1,
                data = data[data$rank_response>10,], family = binomial("logit"))

avg.pred.social.desirability <- predict(fit.list,
                                        direct.glm = fit.sens, se.fit = TRUE)

pred.bias <- as.matrix(avg.pred.social.desirability[[1]])
pred.bias <- data.frame(pred.bias)
pred.bias$est <- c("List \n estimate","Direct \n estimate","Difference")
pred.bias$est <- factor(pred.bias$est)
pred.bias$est = factor(pred.bias$est,levels(pred.bias$est)[c(2,3,1)])

pct10.gr <- ggplot(pred.bias,aes(x=est,y=fit,colour=est)) +
  geom_errorbar(aes(ymin=lwr, ymax=upr), width=.1,position = position_dodge(width = 0.2), size = 0.25) +
  geom_point(position = position_dodge(width = 0.2),size=1.3) +
  scale_y_continuous(limits=c(0,0.5)) +
  labs(y="",x="", title="10th pctl < X") +
  theme_minimal(base_size = 12) + scale_colour_grey(start=0.5, end=0.2) +
  theme(legend.position="none", axis.title.y = element_text(vjust = 2.4),axis.text.x = element_text(size=6.5),
        axis.text.y = element_text(size=5.3),plot.title = element_text(size=9))

fit.list <- ictreg(items ~ 1,
                   data = data[data$rank_response>5&data$rank_response<95,], treat = "list_treatment", J = 4, method = "lm")

fit.sens <- glm(bribe ~ 1,
                data = data[data$rank_response>5&data$rank_response<95,], family = binomial("logit"))

avg.pred.social.desirability <- predict(fit.list,
                                        direct.glm = fit.sens, se.fit = TRUE)

pred.bias <- as.matrix(avg.pred.social.desirability[[1]])
pred.bias <- data.frame(pred.bias)
pred.bias$est <- c("List \n estimate","Direct \n estimate","Difference")
pred.bias$est <- factor(pred.bias$est)
pred.bias$est = factor(pred.bias$est,levels(pred.bias$est)[c(2,3,1)])

pct595.gr <- ggplot(pred.bias,aes(x=est,y=fit,colour=est)) +
  geom_errorbar(aes(ymin=lwr, ymax=upr), width=.1,position = position_dodge(width = 0.2), size = 0.25) +
  geom_point(position = position_dodge(width = 0.2),size=1.3) +
  scale_y_continuous(limits=c(0,0.5)) +
  labs(y="",x="", title="5th pctl < X < 95th pctl") +
  theme_minimal(base_size = 12) + scale_colour_grey(start=0.5, end=0.2) +
  theme(legend.position="none", axis.title.y = element_text(vjust = 2.4),axis.text.x = element_text(size=6.5),
        axis.text.y = element_text(size=5.3),plot.title = element_text(size=9))

# Combine response graphs
list.combined <- ggarrange(pct5.gr, pct10.gr, pct595.gr,
                           ncol = 3, nrow = 1, common.legend = FALSE, legend = "none")
list.combined <- annotate_figure(list.combined, left = text_grob("Estimated proportion",size = 7,rot = 90,hjust = 0.2,vjust = 3))
ggsave("listresponsetime.pdf",width=18,height=6.5, units = "cm")

#####################################################################################################################
# REGRESSION ESTIMATES, APPENDIX
#####################################################################################################################
# TABLE 5
# Regression table, covariates
fit.soc3 <- lm(corr_soc ~ gov_support_01 + corr_soc_treatment + gov_support_01:corr_soc_treatment +
               female + city_over200k + university_ed + income_top20 + Age + I(Age^2), data=data)
fit.econ3 <- lm(econ ~ gov_support_01 + econ_treatment + gov_support_01:econ_treatment +
                  female + city_over200k + university_ed + income_top20 + Age + I(Age^2), data=data)
fit.pol3 <- lm(corr_pol ~ gov_support_01 + corr_pol_treatment + gov_support_01:corr_pol_treatment +
                 female + city_over200k + university_ed + income_top20 + Age + I(Age^2), data=data)
fit.worry3 <- lm(corr_worry ~ gov_support_01 + corr_worry_treatment + gov_support_01:corr_worry_treatment +
                   female + city_over200k + university_ed + income_top20 + Age + I(Age^2), data=data)

names(fit.soc3$coefficients)[names(fit.soc3$coefficients) == "gov_support_01"] <- "Government support"
names(fit.econ3$coefficients)[names(fit.econ3$coefficients) == "gov_support_01"] <- "Government support"
names(fit.pol3$coefficients)[names(fit.pol3$coefficients) == "gov_support_01"] <- "Government support"
names(fit.worry3$coefficients)[names(fit.worry3$coefficients) == "gov_support_01"] <- "Government support"

names(fit.soc3$coefficients)[names(fit.soc3$coefficients) == "corr_soc_treatment"] <- "Prime"
names(fit.econ3$coefficients)[names(fit.econ3$coefficients) == "econ_treatment"] <- "Prime"
names(fit.pol3$coefficients)[names(fit.pol3$coefficients) == "corr_pol_treatment"] <- "Prime"
names(fit.worry3$coefficients)[names(fit.worry3$coefficients) == "corr_worry_treatment"] <- "Prime"

names(fit.soc3$coefficients)[names(fit.soc3$coefficients) == "gov_support_01:corr_soc_treatment"] <- "Gov. support x Prime"
names(fit.econ3$coefficients)[names(fit.econ3$coefficients) == "gov_support_01:econ_treatment"] <- "Gov. support x Prime"
names(fit.pol3$coefficients)[names(fit.pol3$coefficients) == "gov_support_01:corr_pol_treatment"] <- "Gov. support x Prime"
names(fit.worry3$coefficients)[names(fit.worry3$coefficients) == "gov_support_01:corr_worry_treatment"] <- "Gov. support x Prime"

# Tablenotes were added using "threeparttable" in LaTeX
stargazer(fit.soc3,fit.econ3,fit.pol3,fit.worry3, title="The effect of the political prime on perceived corruption, with covariates",
          align=TRUE, dep.var.labels.include = TRUE, 
          dep.var.labels = c("Corruption increase","Worse economy","Corruption in politics","Corruption worry"),
          covariate.labels = c("Government support","Prime","Gov. support x Prime","Female","City inhabitant",
                               "University education","Top 20% income", "Age","Age$^2$"),
          order = c(1,2,9,3,4,5,6,7,8,10), omit.stat=c("LL","ser","f","rsq"), no.space=TRUE, column.sep.width = "-2pt",
          se = starprep(fit.soc3,fit.econ3,fit.pol3,fit.worry3), star.cutoffs = c(0.05,0.01,0.001),
          notes = "All models are estimated using OLS. Robust standard errors in parentheses (HC2). `Prime' is an indicator variable equal to 1 if a respondent was asked the political questions before a specific corruption/economy question, and 0 otherwise. $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001.", digits = 2,
          notes.append = FALSE, notes.align = "l", font.size = "footnotesize")

# TABLE 8
# Regression estimates, response time outliers excluded
fit.soc10 <- lm(corr_soc ~ gov_support_01 + corr_soc_treatment + rank_response + 
                  gov_support_01:rank_response + gov_support_01:corr_soc_treatment + corr_soc_treatment:rank_response +
                  gov_support_01:corr_soc_treatment:rank_response, data=data)
fit.soc11 <- lm(corr_soc ~ gov_support_01 + corr_soc_treatment + gov_support_01:corr_soc_treatment, data=data)
anova(fit.soc11,fit.soc10)

fit.econ10 <- lm(econ ~ gov_support_01 + econ_treatment + rank_response + 
                  gov_support_01:rank_response + gov_support_01:econ_treatment + econ_treatment:rank_response +
                  gov_support_01:econ_treatment:rank_response, data=data)
fit.econ11 <- lm(econ ~ gov_support_01 + econ_treatment + gov_support_01:econ_treatment, data=data)
anova(fit.econ11,fit.econ10)

fit.pol10 <- lm(corr_pol ~ gov_support_01 + corr_pol_treatment + rank_response + 
                   gov_support_01:rank_response + gov_support_01:corr_pol_treatment + corr_pol_treatment:rank_response +
                   gov_support_01:corr_pol_treatment:rank_response, data=data)
fit.pol11 <- lm(corr_pol ~ gov_support_01 + corr_pol_treatment + gov_support_01:corr_pol_treatment, data=data)
anova(fit.pol11,fit.pol10)

fit.worry10 <- lm(corr_worry ~ gov_support_01 + corr_worry_treatment + rank_response + 
                  gov_support_01:rank_response + gov_support_01:corr_worry_treatment + corr_worry_treatment:rank_response +
                  gov_support_01:corr_worry_treatment:rank_response, data=data)
fit.worry11 <- lm(corr_worry ~ gov_support_01 + corr_worry_treatment + gov_support_01:corr_worry_treatment, data=data)
anova(fit.worry11,fit.worry10)

fit.soc12 <- lm_robust(corr_soc ~ gov_support_01 + corr_soc_treatment + gov_support_01:corr_soc_treatment, data=data[data$rank_response>5,])
fit.soc13 <- lm_robust(corr_soc ~ gov_support_01 + corr_soc_treatment + gov_support_01:corr_soc_treatment, data=data[data$rank_response>10,])
fit.soc14 <- lm_robust(corr_soc ~ gov_support_01 + corr_soc_treatment + gov_support_01:corr_soc_treatment, data=data[data$rank_response>5&data$rank_response<95,])

fit.econ12 <- lm_robust(econ ~ gov_support_01 + econ_treatment + gov_support_01:econ_treatment, data=data[data$rank_response>5,])
fit.econ13 <- lm_robust(econ ~ gov_support_01 + econ_treatment + gov_support_01:econ_treatment, data=data[data$rank_response>10,])
fit.econ14 <- lm_robust(econ ~ gov_support_01 + econ_treatment + gov_support_01:econ_treatment, data=data[data$rank_response>5&data$rank_response<95,])

fit.pol12 <- lm_robust(corr_pol ~ gov_support_01 + corr_pol_treatment + gov_support_01:corr_pol_treatment, data=data[data$rank_response>5,])
fit.pol13 <- lm_robust(corr_pol ~ gov_support_01 + corr_pol_treatment + gov_support_01:corr_pol_treatment, data=data[data$rank_response>10,])
fit.pol14 <- lm_robust(corr_pol ~ gov_support_01 + corr_pol_treatment + gov_support_01:corr_pol_treatment, data=data[data$rank_response>5&data$rank_response<95,])

fit.worry12 <- lm_robust(corr_worry ~ gov_support_01 + corr_worry_treatment + gov_support_01:corr_worry_treatment, data=data[data$rank_response>5,])
fit.worry13 <- lm_robust(corr_worry ~ gov_support_01 + corr_worry_treatment + gov_support_01:corr_worry_treatment, data=data[data$rank_response>10,])
fit.worry14 <- lm_robust(corr_worry ~ gov_support_01 + corr_worry_treatment + gov_support_01:corr_worry_treatment, 
                         data=data[data$rank_response>5&data$rank_response<95,])

#Table
resp.table <- matrix(nrow = 16,ncol = 3)
resp.table[1,] <- c("","$5th pctl < X$","$10th pctl < X$","$5th pctl < X < 95th$")
resp.table[,1] <- c("","Corr. increase $\\beta_{1}$","","Corr. increase $\\delta$","","Worse economy $\\beta_{1}$","","Worse economy $\\delta$","",
                    "Corr. in politics $\\beta_{1}$","","Corr. in politics $\\delta$","","Corr. worry $\\beta_{1}$","","Corr. worry $\\delta$","")

resp.table[1,1] <- round(fit.soc12$coefficients[[2]],3)
resp.table[3,1] <- round(fit.soc12$coefficients[[4]],3)
resp.table[2,1] <- round(fit.soc12$std.error[[2]],3)
resp.table[4,1] <- round(fit.soc12$std.error[[4]],3)

resp.table[1,2] <- round(fit.soc13$coefficients[[2]],3)
resp.table[3,2] <- round(fit.soc13$coefficients[[4]],3)
resp.table[2,2] <- round(fit.soc13$std.error[[2]],3)
resp.table[4,2] <- round(fit.soc13$std.error[[4]],3)

resp.table[1,3] <- round(fit.soc14$coefficients[[2]],3)
resp.table[3,3] <- round(fit.soc14$coefficients[[4]],3)
resp.table[2,3] <- round(fit.soc14$std.error[[2]],3)
resp.table[4,3] <- round(fit.soc14$std.error[[4]],3)

resp.table[5,1] <- round(fit.econ12$coefficients[[2]],3)
resp.table[7,1] <- round(fit.econ12$coefficients[[4]],3)
resp.table[6,1] <- round(fit.econ12$std.error[[2]],3)
resp.table[8,1] <- round(fit.econ12$std.error[[4]],3)

resp.table[5,2] <- round(fit.econ13$coefficients[[2]],3)
resp.table[7,2] <- round(fit.econ13$coefficients[[4]],3)
resp.table[6,2] <- round(fit.econ13$std.error[[2]],3)
resp.table[8,2] <- round(fit.econ13$std.error[[4]],3)

resp.table[5,3] <- round(fit.econ14$coefficients[[2]],3)
resp.table[7,3] <- round(fit.econ14$coefficients[[4]],3)
resp.table[6,3] <- round(fit.econ14$std.error[[2]],3)
resp.table[8,3] <- round(fit.econ14$std.error[[4]],3)

resp.table[9,1] <- round(fit.pol12$coefficients[[2]],3)
resp.table[11,1] <- round(fit.pol12$coefficients[[4]],3)
resp.table[10,1] <- round(fit.pol12$std.error[[2]],3)
resp.table[12,1] <- round(fit.pol12$std.error[[4]],3)

resp.table[9,2] <- round(fit.pol13$coefficients[[2]],3)
resp.table[11,2] <- round(fit.pol13$coefficients[[4]],3)
resp.table[10,2] <- round(fit.pol13$std.error[[2]],3)
resp.table[12,2] <- round(fit.pol13$std.error[[4]],3)

resp.table[9,3] <- round(fit.pol14$coefficients[[2]],3)
resp.table[11,3] <- round(fit.pol14$coefficients[[4]],3)
resp.table[10,3] <- round(fit.pol14$std.error[[2]],3)
resp.table[12,3] <- round(fit.pol14$std.error[[4]],3)

resp.table[13,1] <- round(fit.worry12$coefficients[[2]],3)
resp.table[15,1] <- round(fit.worry12$coefficients[[4]],3)
resp.table[14,1] <- round(fit.worry12$std.error[[2]],3)
resp.table[16,1] <- round(fit.worry12$std.error[[4]],3)

resp.table[13,2] <- round(fit.worry13$coefficients[[2]],3)
resp.table[15,2] <- round(fit.worry13$coefficients[[4]],3)
resp.table[14,2] <- round(fit.worry13$std.error[[2]],3)
resp.table[16,2] <- round(fit.worry13$std.error[[4]],3)

resp.table[13,3] <- round(fit.worry14$coefficients[[2]],3)
resp.table[15,3] <- round(fit.worry14$coefficients[[4]],3)
resp.table[14,3] <- round(fit.worry14$std.error[[2]],3)
resp.table[16,3] <- round(fit.worry14$std.error[[4]],3)

row.names(resp.table) <- c("Corr. increase beta_{1}","","Corr. increase delta","","Worse economy beta_{1}","","Worse economy delta","",
                           "Corr. in politics beta_{1}","","Corr. in politics delta","","Corr. worry beta_{1}","","Corr. worry delta","")

colnames(resp.table) <- c("5th pctl < X","10th pctl < X","5th pctl < X < 95th")

for (i in c(2,4,6,8,10,12,14,16)) {
  for(j in 1:3){
    resp.table[i,j] <- paste("(",resp.table[i,j],")", sep="")
  }
}

# Tablenotes were added using "threeparttable" in LaTeX
stargazer(resp.table, 
          notes = "\\beta_{1} and \\delta refer to coefficient estimates in equation (1). All models estimated using OLS. X denotes response time percentiles. Robust standard errors in parentheses (HC2).",
          title = "Estimates of eq. (1), excluding response-time outliers", column.sep.width = "-2pt",font.size = "footnotesize", digits = 2)

# TABLE 7
# Regression table, gov support (3 cat)
data$gov_supp_3cat <- factor(data$gov_supp_3cat,labels = c("Opposing gov.","Neither","Favoring gov."))
fit.soc4 <- lm(corr_soc ~ gov_supp_3cat + corr_soc_treatment + gov_supp_3cat:corr_soc_treatment, data=data)
fit.econ4 <- lm(econ ~ gov_supp_3cat + econ_treatment + gov_supp_3cat:econ_treatment, data=data)
fit.pol4 <- lm(corr_pol ~ gov_supp_3cat + corr_pol_treatment + gov_supp_3cat:corr_pol_treatment, data=data)
fit.worry4 <- lm(corr_worry ~ gov_supp_3cat + corr_worry_treatment + gov_supp_3cat:corr_worry_treatment, data=data)

names(fit.soc4$coefficients)[names(fit.soc4$coefficients) == "gov_supp_3catNeither"] <- "Government: Neutral"
names(fit.econ4$coefficients)[names(fit.econ4$coefficients) == "gov_supp_3catNeither"] <- "Government: Neutral"
names(fit.pol4$coefficients)[names(fit.pol4$coefficients) == "gov_supp_3catNeither"] <- "Government: Neutral"
names(fit.worry4$coefficients)[names(fit.worry4$coefficients) == "gov_supp_3catNeither"] <- "Government: Neutral"

names(fit.soc4$coefficients)[names(fit.soc4$coefficients) == "gov_supp_3catFavoring gov."] <- "Government: Favor"
names(fit.econ4$coefficients)[names(fit.econ4$coefficients) == "gov_supp_3catFavoring gov."] <- "Government: Favor"
names(fit.pol4$coefficients)[names(fit.pol4$coefficients) == "gov_supp_3catFavoring gov."] <- "Government: Favor"
names(fit.worry4$coefficients)[names(fit.worry4$coefficients) == "gov_supp_3catFavoring gov."] <- "Government: Favor"

names(fit.soc4$coefficients)[names(fit.soc4$coefficients) == "corr_soc_treatment"] <- "Prime (Government: Oppose)"
names(fit.econ4$coefficients)[names(fit.econ4$coefficients) == "econ_treatment"] <- "Prime (Government: Oppose)"
names(fit.pol4$coefficients)[names(fit.pol4$coefficients) == "corr_pol_treatment"] <- "Prime (Government: Oppose)"
names(fit.worry4$coefficients)[names(fit.worry4$coefficients) == "corr_worry_treatment"] <- "Prime (Government: Oppose)"

names(fit.soc4$coefficients)[names(fit.soc4$coefficients) == "gov_supp_3catNeither:corr_soc_treatment"] <- "Government: Neutral x Prime"
names(fit.econ4$coefficients)[names(fit.econ4$coefficients) == "gov_supp_3catNeither:econ_treatment"] <- "Government: Neutral x Prime"
names(fit.pol4$coefficients)[names(fit.pol4$coefficients) == "gov_supp_3catNeither:corr_pol_treatment"] <- "Government: Neutral x Prime"
names(fit.worry4$coefficients)[names(fit.worry4$coefficients) == "gov_supp_3catNeither:corr_worry_treatment"] <- "Government: Neutral x Prime"

names(fit.soc4$coefficients)[names(fit.soc4$coefficients) == "gov_supp_3catFavoring gov.:corr_soc_treatment"] <- "Government: Favor x Prime"
names(fit.econ4$coefficients)[names(fit.econ4$coefficients) == "gov_supp_3catFavoring gov.:econ_treatment"] <- "Government: Favor x Prime"
names(fit.pol4$coefficients)[names(fit.pol4$coefficients) == "gov_supp_3catFavoring gov.:corr_pol_treatment"] <- "Government: Favor x Prime"
names(fit.worry4$coefficients)[names(fit.worry4$coefficients) == "gov_supp_3catFavoring gov.:corr_worry_treatment"] <- "Government: Favor x Prime"

# Tablenotes were added using "threeparttable" in LaTeX
stargazer(fit.soc4,fit.econ4,fit.pol4,fit.worry4, title="The effect of the political prime on perceived corruption, oppose/favor government",
          align=TRUE, dep.var.labels.include = TRUE, 
          dep.var.labels = c("Corruption increase","Worse economy","Corruption in politics","Corruption worry"),
          omit.stat=c("LL","ser","f","rsq"), no.space=TRUE, se = starprep(fit.soc4,fit.econ4,fit.pol4,fit.worry4),
          notes = "All models are estimated using OLS. Robust standard errors in parentheses (HC2). `Prime' is an indicator variable equal to 1 if a respondent was asked the political questions before a specific corruption/economy question, and 0 otherwise. `Opposing' corresponds to category 1-2 on the government support variable, `neutral' corresponds to category 3, and `favoring' corresponds to category 4-5. $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001.", 
          notes.append = FALSE, notes.align = "l", column.sep.width = "-1pt", digits = 2,
          order = c(3,1,2,4,5), star.cutoffs = c(0.05,0.01,0.001), font.size = "footnotesize")

# TABLE 6
# Regression table, OLR
data2 <- data

data2$corr_soc <- factor(data2$corr_soc,labels = c("Decreased a lot","Decreased somewhat","Stayed the same","Increased somewhat","Increased a lot"))
data2$econ <- factor(data2$econ,labels = c("Got a lot stronger","Got a little strongert","Stayed the same","Got a little weaker","Got a lot weaker"))
data2$corr_pol <- factor(data2$corr_pol,labels = c("Almost none","A few","Some","Quite a lot","Almost all"))
data2$corr_worry <- factor(data2$corr_worry,labels = c("Not worried at all","A little worried","Somewhat worried","Very worried"))

fit.soc6 <- polr(corr_soc ~ gov_support_01 + corr_soc_treatment + gov_support_01:corr_soc_treatment, data=data2, Hess = TRUE)
fit.econ6 <- polr(econ ~ gov_support_01 + econ_treatment + gov_support_01:econ_treatment, data=data2, Hess = TRUE)
fit.pol6 <- polr(corr_pol ~ gov_support_01 + corr_pol_treatment + gov_support_01:corr_pol_treatment, data=data2, Hess = TRUE)
fit.worry6 <- polr(corr_worry ~ gov_support_01 + corr_worry_treatment + gov_support_01:corr_worry_treatment, data=data2, Hess = TRUE)

names(fit.soc6$coefficients)[names(fit.soc6$coefficients) == "gov_support_01"] <- "Government support"
names(fit.econ6$coefficients)[names(fit.econ6$coefficients) == "gov_support_01"] <- "Government support"
names(fit.pol6$coefficients)[names(fit.pol6$coefficients) == "gov_support_01"] <- "Government support"
names(fit.worry6$coefficients)[names(fit.worry6$coefficients) == "gov_support_01"] <- "Government support"

names(fit.soc6$coefficients)[names(fit.soc6$coefficients) == "corr_soc_treatment"] <- "Prime"
names(fit.econ6$coefficients)[names(fit.econ6$coefficients) == "econ_treatment"] <- "Prime"
names(fit.pol6$coefficients)[names(fit.pol6$coefficients) == "corr_pol_treatment"] <- "Prime"
names(fit.worry6$coefficients)[names(fit.worry6$coefficients) == "corr_worry_treatment"] <- "Prime"

names(fit.soc6$coefficients)[names(fit.soc6$coefficients) == "gov_support_01:corr_soc_treatment"] <- "Gov. support x Prime"
names(fit.econ6$coefficients)[names(fit.econ6$coefficients) == "gov_support_01:econ_treatment"] <- "Gov. support x Prime"
names(fit.pol6$coefficients)[names(fit.pol6$coefficients) == "gov_support_01:corr_pol_treatment"] <- "Gov. support x Prime"
names(fit.worry6$coefficients)[names(fit.worry6$coefficients) == "gov_support_01:corr_worry_treatment"] <- "Gov. support x Prime"

# Tablenotes were added using "threeparttable" in LaTeX
stargazer(fit.soc6,fit.econ6,fit.pol6,fit.worry6, title="The effect of the political prime on perceived corruption, OLR estimates",
          align=TRUE, dep.var.labels.include = TRUE, 
          dep.var.labels = c("Corruption increase","Worse economy","Corruption in politics","Corruption worry"),
          omit.stat=c("LL"), no.space=TRUE, model.names = FALSE,
          notes = "All models are estimated using ordinal logistic regression. Standard errors in parentheses. `Prime' is an indicator variable equal to 1 if a respondent was asked the political questions before a specific corruption/economy question, and 0 otherwise. $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001.",  
          notes.append = FALSE, notes.align = "l", column.sep.width = "-2pt", digits = 2,
          star.cutoffs = c(0.05,0.01,0.001), font.size = "footnotesize")
